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This paper describes a modeling method and a new optimal control approach to investigate a Mach number 
control problem for the NASA Ames 11-Foot Transonic Wind Tunnel. The flow in the wind tunnel is modeled 
by the 1-D unsteady Euler equations whose boundary conditions prescribe a controlling action by a compres- 
sor. The boundary control inputs to the compressor are in turn controlled by a drive motor system and an inlet 
guide vane system whose dynamics are modeled by ordinary differential equations. The resulting Euler equa- 
tions are thus coupled to the ordinary differential equations via the boundary conditions. Optimality conditions 
are established by an adjoint method and are used to develop a model predictive linear-quadratic optimal con- 
trol for regulating the Mach number due to a test model disturbance during a continuous pitch. 


Introduction 

In this paper, we present a new optimal control approach for a continuous system governed by first order hyperbolic 
partial differential equations (PDEs) coupled to a discrete system modeled by ordinary differential equations (ODEs) 
imposed at the continuous system boundary. Hyperbolic equations model many physical systems which are governed 
by conservation laws such as fluid flow, and therefore have been thoroughly studied. However, problems involving 
hyperbolic PDEs coupled to ODEs via boundary conditions are much less common. The goal of this paper is to 
present a model predictive optimal control for solving a Mach number control problem for a wind tunnel using this 
new modeling approach. 

Fig. 1 illustrates the NASA Ames 1 1-Ft Transonic Wind Tunnel. It is capable of delivering an air speed from Mach 
0.2 to Mach 1 .5 in a test section measured 1 1 ft wide by 1 1 ft high. The air flow is delivered by a compressor driven by 
a set of synchronous induction AC motors operated in a speed range from 310 to 645 rpm with a total input power of 
236,000 horsepowers. The compressor is a variable-geometry machine equipped with inlet guide vanes (IG Vs) having 
adjustable trailing edge flaps through a range of angular deflections from —7.5° to 19.5°. 



Fig. 1 - NASA Ames 11 -Foot Transonic Wind Tunnel 
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The test section aerodynamic condition is normally controlled to ensure that the Mach number variation is within 
a prescribed tolerance, which can be as small as 0.001. For subsonic flow, the Mach number control is accomplished 
by adjusting the IGV flap deflection while holding the compressor speed constant at the discrete set points as shown 
in Fig. 2 by a compressor speed feedback control. 



Fig. 2 - NASA Ames 1 1-Ft Transonic Wind Tunnel Mach Envelope 

Typically, during a wind tunnel test, the test model is pitched through a series of angles of attack at a nominal 
Mach number. This is known as a pitch polar. As the angle of attack increases, the drag force on the test model also 
increases. This drag force thus translates into a loss in the total pressure across the test model. The total pressure 
disturbance then propagates downstream to the compressor inlet and sets up a dynamic imbalance in the compressor 
equilibrium, thereby causing the compressor speed to fall. This consequently results in a deviation in the Mach number 
in the test section from the desired set point. Presently, to minimize this Mach number deviation, the test model has 
to be paused in between changes in the pitch angle so that the flow condition can be reestablished by a Mach number 
feedback control using the IGVs. 

It is a recognized that the current Mach number feedback control strategy is deficient and therefore an alternate 
Mach number control strategy should be considered to provide an effective Mach number control strategy that will 
allow the test model to be pitched continuously. This strategy would translate into a significant advantage over the 
current pitch-pause mode. Therefore, the motivation of this study is to investigate a model predictive Mach number 
optimal control for a wind tunnel based on a modeling approach of coupled hyperbolic PDEs and ODEs. 

Wind Tunnel Modeling 

Wind tunnel flow modeling for control design purposes using unsteady fluid dynamics equations has rarely been 
considered because of the inherent mathematical complexity. Understandably, many wind tunnel control models have 
been based on heuristic or ad hoc methods. Soeterboek describes a method for control modeling of a test section 
Mach number by an experimental transfer function coupled with a time delay 1 . The current Mach number control 
for NASA Ames 1 1-Ft Transonic Wind Tunnel was developed using a quasi- steady state model based on an empirical 
relationship between the Mach number and the compressor performance characteristics 2 . The quasi-steady state model 
circiimvents_the need for modeling the complex unsteady _air flow, but does not capture the time delay effect. As a 
result, it can lead to poor control handling of disturbance rejection due to the fluid transport delay that exists between 
the test section and the compressor in a typical wind tunnel. 

In this paper, we propose to establish a model of the wind tunnel flow based on the governing conservation laws of 
physics to capture the fluid transiency. In the process, this model should be more accurate than a heuristic model and 

• 

2 

American Institute of Aeronautics and Astronautics 




thus can be used for the proposed predictive Mach number optimal control strategy. The governing Euler equations 
for an unsteady 1-D flow are 3 

|(M) + J^) = ° (i) 

I {puA) + 1 ( pA + ^ - pd i + l pu 2 A i = 0 (2) 

~ ^peA + ^pu 2 A^j + ~ ^ pueA + ^pu 3 A + puAxj = 0 (3) 

where pis the static density, u is the flow speed, p is the static pressure, e is the internal energy, A is the cross-sectional 
area, D is the hydraulic diameter, and / is the friction factor. 

Equations (l)-(3) are the conservation equations for fluid flow. It is, however, usually more convenient to express 
these equations in terms of the mass flow m, the total pressure po, and the total temperature To. The Mach number 
can be related to these three quantities by 

i fe+i 

™ = \fw 0 PoAM ( x + m2 ) ^ (4) 

where M is the Mach number and k is the specific heat ratio. 

We now introduce an alternate form of the Euler equations as follows 

y t + A (y, x) y x + B (y , x) = 0 (5) 

, . r ~iT 

where y (x, t) = [ m p 0 To ] ~ is a flow vector and 
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where c is the speed of sound, T is the static temperature, and po is the total density. 

The flow recirculation in a wind tunnel is controlled by a compressor, which creates a pressure rise to compensate 
for the pressure losses throughout the wind tunnel circuit. We let x = 0 and x = L be the compressor exit and inlet 
stations, respectively. Then the boundary conditions for Eq. (5) are specified by the following compressor performance 
model, which can be obtained v empirically from experimental data as 

m (0, t) = m ( L , t) (6) 
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(7) 

( 8 ) 


where bi , <%, are the empirical coefficients. 

Equations (6)-(8) are the respective continuity, momentum, and energy equations that relate the mass flow, the total 
pressure rise, and the temperature rise across the compressor to the compressor speed a; and the IGV flap deflection 0 
which are the control inputs to the compressor. These equations can be written generally as 


y (0, t) = g (y {L, t) , u) (9) 

where u = [ 0 (J ] T is called a boundary control vector. 

The wind tunnel usually operates from one steady state condition to another. Therefore, the initial condition for 
Eq. (5) is specified as follows 


y(x,0) = h(x) 
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where h (x) is a steady state solution of Eq. (5), 

Since the compressor speed must be controlled simultaneously with the Mach number, we must take into account 
the dynamics of the drive motors and the IGV system that actually control the compressor speed and IGV flap de- 
flection. The three-phase synchronous induction AC drive motors use a liquid rheostat system as a primary means of 
changing the rotor resistance R r which changes the motor torque to control the compressor speed. The dynamics of 
this drive motor system is described by the following motor torque equation 


oj ■■ 


KqjiRpUJs {^S ^0 


[R s {oJ s - u) + Rr^s] + (u) s ~ W) 


■ K a [po (0, t ) - Po ( L , i)] 


(ID 


where R s is the stator resistance, L s is the stator inductance, oj s is the synchronous speed, and K m and K a are the 
torque constants. 

The IGVs are driven by DC servo motors and controlled by the motor field voltage V a . The dynamics of this 
system is modeled as 

e = K v Va - \p (L, t) u 2 (L, t) (Ke6 + K c ) (12) 

where K v , Kg, K c are the IGV system constants. 

We observe that the boundary control inputs oj and 6 to the compressor become the states or outputs of the drive 
motors and the IGV system. Both Eqs. (11) and (12) are essentially the actuator dynamics of the compressor and are 
observed to also depend on the flow variables at the boundary. Thus, these equations are coupled with Eq. (5) through 
the boundary conditions (6) to (8) and can be expressed in a general form as 


U = f (y (0, t) , y (L, t) , u, v) 


(13) 


where v = [ V a R r ] T is a control vector. 

For a well-posed problem, the initial condition for Eq. (13) must be consistent with the initial condition (10) and 
the boundary condition (9) such that 

h (0) = g (h (L) , u (0)) (14) 


Computationai Method 


It can be shown that the eigenvalues of the matrix A are A (A) = u,u ± c, which are the wave propagation 
speeds in a fluid medium. Subsonic flow therefore involves two waves propagating downstream and one wave prop- 
agating upstream from the source. If A + = diag ( u 4- c> u, 0) is a diagonal matrix of the positive eigenvalues, 
A“ = diag (0, 0, u — c) is a diagonal matrix of the negative eigenvalues, and $ is a matrix of the right eigenvec- 
tors of the matrix A, then the matrix A can be split into a semi-positive definite matrix and a semi-negative definite 
matrix as follows 

A — A + -h A“ (15) 


where 


A+ = $A + #“ a 
A” = ^A”#” 1 


Equation (5) can he cast into a discretized characteristic form using a combined first order upwind and downwind 
explicit finite-difference method 4 as follows 


A + $ _1 A t 

*4^+1 = y« - (y a - y *-w) - ^ (y wj - y«) - 


(16) 


where i = 2, 3, . . . , m — 1 denotes the interior points not on the boundary, j = 2, 3, . . . , n — 1 denotes the time step. 
Ax and At. are the space and time steps which must satisfy the following CEL condition 4 


At < 


Ax 
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( 17 ) 
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For subsonic flow, we require the number of boundary conditions at x = 0 be equal to the number of positive 
eigenvalues, or two. We then combine the two boundary conditions (7) and (8) with the characteristic equation in the 
third row of Eq. (16) to compute the flow variables at x = 0 as follows 
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where I* is the z-th row of a 3-by-3 identity matrix, 0i X 3 is a l-by-3 zero vector, and = #“ 1 . 

At x = L, we combine the remaining boundary condition (6) with the two remaining characteristic equations from 
Eq. (16) as follows 
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(19) 


Both Eqs. (18) and (19) are nonlinear and are coupled together. An iterative method is used to solve for the 
flow conditions at the compressor inlet and exit simultaneously. Once the solution at the boundary is established, the 
solution for the interior points can be computed in a time marching procedure. 


Validation 

To validate the wind tunnel distributed model, we compute a transition from Mach 0.6 at 455 rpm to Mach 0.9 at 
590 rpm. The wind tunnel distributed model is discretized into 45 nodes as shown in Fig. 3. 



0 100 200 300 400 500 600 700 

x, ft 


Fig. 3 - Wind Tunnel Discretization 

The step size in x is chosen so that important features of the wind tunnel such as the aftercooler and the test section 
dimensions are captured in the discretization. This step size is found to be Ax — 18.07 ft. The time step is then 
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chosen to be At — 0.01 sec in order to satisfy the CFL condition (17). At each node, there are three flow quantities, so 
the discretized model contains a total of 135 aerodynamic states plus two compressor dynamic states and two control 
inputs. 

The time histories of the compressor speed and the IGV flap deflection during the Mach number transition are 
shown in Fig. 4. The IGV flap deflection is first adjusted to its maximum value while the compressor speed is 
regulated at 455 rpm. A new compressor speed set point command at 590 rpm is then sought while the IGV flap 
deflection is maintained at its maximum value. When the compressor settles to the new set point, the IGV flaps are 
then actuated until the Mach number reaches 0.9. A small ripple in the compressor speed occurs at the beginning of 
the final IGV flap actuation due to an integral feedback control regulating the motor speed set point as a result of the 
changing aerodynamic torque during the IGV actuation. 

Fig. 5 shows the resulting Mach number variation throughout the transition. The computed Mach number from 
the distributed model agrees very well with the experimental data, thus validating the modeling approach. 



Fig. 4 - Compressor Control Time History 



Fig. 5 - Mach Number Response 
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Fig. 6 is the plot of the Mach number distribution along the wind tunnel showing the Mach number response in 
the test section. The demonstrated accuracy of the wind tunnel model is critical to the model predictive Mach number 
control strategy wherein the control inputs are strictly computed from the mathematical model of the wind tunnel. 



Fig. 6 - Mach Number Distribution 


Adjoint Method 

In order to develop a model predictive Mach number optimal control, we propose to establish an optimal control 
theory for systems modeled by PDEs coupled to ODEs via boundary conditions based on an adjoint method 5 . Optimal 
control for PDE systems has been examined in the literature 6,7,8 , but little existing theory has been found for the 
systems of the present interest. The adjoint method has been used in many fluid optimization problems 9,10 . There are 
two approaches to the adjoint formulation in the optimization application domain. This first approach is to formulate 
the adjoint method directly to the fluid equations of motion. This is called a continuous adjoint method 11 . The 
other approach is to formulate the adjoint method to the discretized equations of motion and as such is called a 
discrete adjoint method. Studies have found that the continuous adjoint method is generally preferred 11 . We adopt 
the continuous adjoint approach in this study. Towards that end, we are interested in obtaining an optimal control that 
will minimize the responses at the system boundaries and the boundary control action by considering the following 
linear-quadratic cost function with a fixed terminal time t / 

min J(y(0,t),y(L,t),u,v) - [y T (0, t) Py (0, t) + y T (L, t) Sy (L, t) + u T Qu + v T Rv] dt (20) 

subject to Eq. (5) coupled to Eq. (12) via the boundary condition (9) , where P, S and Q are semi-positive definite 
symmetric matrices, and R is a positive definite symmetric matrix. 

We define A (x, t) as the adjoint vector for Eq. (5) and fi (£) as the adjoint vector for Eq. (12). Using the Lagrange 
multiplier method, the augmented cost function becomes 


J (y (0, t) , y (L, t ) , u, v) = - / f X T (y t + Ay x + B) 

Jo Jo 

+ f |^y T (o, t) Py (0, t) + iy T (L, t) Sy (L, t) + iu T Qu + iv T Rv + M (~u + f) 


) dxdt 


dt (21) 
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To establish the necessary conditions of optimality, we compute the first variation of the augmented cost function 
using the Green’s theorem which results in 

SJ = [ f 5y T [A t 4- (A T X) x - By \] T dxdt + f 5y T (x,tf) A (x,tf) dx 
Jo J o Jo 

- / [<*y T (L, t ) A T (y (L,t) , L) A (L, t) - dy T (0, t ) A T (y (0, t) , 0) A (0, i)] dt 

J o 

+ (<5y T (0, t) [Py (0, t) + fj (0>t) £t] + 5y T ( L , t) [Sy (L, t) + 

+<Ju T (Qu + £ TfJ, + A) + (Rv + } dt - 5 u t ( tf ) {i ( tf ) (22) 

In addition, we also compute the first variation of the boundary condition (9) as follows 

5y T (0, t) = Sy T (L, t) Sy (Lt t) + (23) 

Substituting Eq. (23) into Eq. (22) and setting 5 J = 0, we obtain the following adjoint PDE 

A) + (A t A) x - By A = 0 (24) 

subject to a boundary condition 

A t (y (L, t),L) A (L, t ) = gy (£ , t) A T (y (0, i) , 0) A (0,i)+ + gy( L , t) f^ 0l t)] M+gya,t)Py (0, t)+ Sy (£, 4) 

(25) 

and a terminal time condition 

A(*,t / ) = 0 (26) 

In addition, we obtain the following adjoint ODE 

A = — gu A T (y (0,t) , 0) A (0, t) - [fj + gufy( 0 ,t)] M - Qu - g^Py (0, i) (27) 

subject to a terminal time condition 

U(t f ) - 0 (28) 

Equations (24) to (28) constitute an adjoint system for the original problem defined by Eqs. (5), (9), (10), (12), 
and (14). These systems represent a duality wherein the adjoint variables are mapped into space and time reversal 
coordinates. Solutions for such systems are generally complex for nonlinear problems and are usually solved by 
gradient methods. For linear problems, they can be solved by other methods which will be demonstrated by applying 
the adjoint method to the Mach number control problem. 

Model Predictive Mach Number Optimal Control 


The optimal control results can now be applied to design a predictive Mach number control for a wind tunnel. One 
of the objectives of a wind tunnel Mach number control is to control the Mach number in the test section due to a 
small total pressure disturbance created by the test model drag as accurately as possible. Neglecting the continuity and 
energy equations since the mass flow and total temperature are less sensitive to the total pressure disturbance caused 
by the test model than the total pressure, then linearizing Eq. (5) about a steady state nominal operating condition 
results in a linearized momentum equation 

Sjbft + <»> 


where y (, x , t) — A p 0 (x, t) is the total pressure perturbation and 


a(x) = M ( x ) c (x) 


4 — 2fc + (fc — 1) M 2 (x) 
2 4- (& — 1) M 2 (x) 
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P (x, t ) = 


kM 2 {x)[l + kM 2 {x)] 

2 [l — M 2 (a:)] _ j 

vj (t) s= ~p 0 ( 3 ; t) M 2 (»t) 


/(g) , gp (t) An 
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where the overbar denotes the nominal steady state condition and w (t) is the total pressure disturbance as a function 
of the nominal Mach number M in the test section at x — xt, the test model drag coefficient Cp, the test model 
reference area A m , and the reference length L m . The validity of this equation excludes the transonic region when 
M {xt) =s 1. 

Linearizing Eq. (7) results in the following boundary condition 


y (0, t) = G Au (t) 4- Hy (L, t) 


(30) 


QT» 

where Au = [ AO Auj £ ] is a compressor boundary control error vector comprising the IGV flap deflection 
error AO, the compressor speed error Au/, and the compressor speed error integral £ = Aujdr. The error integral is 
designed to achieve a zero steady state error in the compressor speed. In essence, the control scheme is a proportional- 

integral control. The matrix G = [ dpo d ^ 0 j and H — are evaluated from Eq. (7). 

The compressor boundary control error dynamics is obtained by linearizing Eq. (13) 


Ah = C Au + D Av + E y (0, t) + F y (L, t) (3 1) 

where Av = [ AV a A R r ] T is the control perturbation vector and 
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To design a Mach number control, we recognize that the Mach number perturbation is caused by a total pressure 
disturbance downstream of the test model. Since there is no total pressure disturbance upstream of the test model, to 
minimize the Mach number perturbation, we would need to minimize the total pressure perturbation at the compressor 
exit at x = 0. Thus, we consider the following linear-quadratic cost function 


J 


-f 


\py 2 (0, t) + iAu r QAu + Av t RAv 
z z z 


dt 


where P > 0, Q > 0, and R > 0. 

Applying the necessary conditions of optimality results in the following adjoint equations 


dX d t vx ^ x 
dt + dx iaX) ~ pa A_0 

a (L) A (L, t ) = Ha (0) A (0, f) + HPy (0, f) + (F t + HE T ) n 
ii = -QAu - (C T + G 7 E t ) fj, - G T a (0) A (0, t) - G T Py (0, t) 
RAv + D r /it = 0 

It can be shown that the general solution of Eq. (29) is 


y (x, f) - a (x,t) [f (f - t d (x)) - q ( x , f)] 


(32) 


(33) 

(34) 

(35) 

(36) 

(37) 


where 



da 
a {a) 


a (X, t) = e~ So ^,t~td(x)+t d (a))da 
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q (x, t) = l a 1 (a, t- td ( x ) + td (<?)) V) (t - td (x) + td (a)) da 

J o 

We now consider two special eases: a short time horizon when tf < tjo , where tr, = td ( L ) is a fluid transport 
delay, and a long time horizon when tf — * oo. For a short time horizon, we have the following relationship 


y (0, t) = GAu (f) - /la (I, f) g (I, t) (38) 

V ( L > t) ~ -a (L, t) q (L, t) (39) 

a (0) A (0, t) — 0 (40) 

On the other hand, for a long time horizon, we obtain the following quasi-steady state solution 

2/(0, t) — [1 — Ha (L, t )] -1 [GAu (t) — Ha (L, t) q (L, f)] (41) 

y (L, t) = a (If, t) [1 - Ha ( L , t)j _1 [GAu (t) - Ha (I, t) q (L, t)] * a (I, f) g (L, t) (42) 

a (0) A (0, t) = a (I, t) [1 - Ha (L, t )] -1 [. HPy (0, t ) + (F t + i?E T ) /a (f)] (43) 

We now define a time horizon parameter 7 which takes on a value ranging from 0 to 1 that represents the effect 
of the time horizon on the optimal control solution such that 7 = 0 corresponds to a short time horizon and 7 = 1 
corresponds to a long time horizon. The optimal control solution can now be obtained from the following the Riccati 
equation by letting n = WAu + V 

W + WC e + c Jw - WDR- 1 D r W + Qe = 0 (44) 

V + C e T V - WDR _1 D t V - WS - U = 0 (45) 


where 

C e = C + EG + 7 (F + EH) Ga (I, i) [1 - jHa (I, t )] _1 
Q e = Q + G t PG [1 - jHa (I, f)] " 2 
S = (F + EH) a (L, t) q (L, t) [1 - ~fHa (L, t)\~ l 
U = G T PHa ( L , t) q (I, t) [l - 7 Ha (I, t )]~ 2 
The control perturbation Av therefore is expressed as a feedback control as 

Av (f) = -R -1 D r W ( t ) Au (t) - R" x D t V ( t ) (46) 


Equation (46) compute the corrective control inputs for the drive motors and the IGV system which are to be added 
to the steady state control inputs at a nominal Mach number to give the total control inputs for the drive motors and 
IGV system according to 


V«{t) 
R r ( t ) 



A V a (t) ' 
A & (t) _ 


(47) 


However, we recognize that the compressor generally does not sense the effect of the total pressure disturbance 
signal caused by the motion of the test model at the same time instance as the disturbance itself. Because of the fluid 
transport delay, there is a time lag between the test section and the compressor inlet such that the compressor dynamic 
response lags the total pressure disturbance in the test section by a time delay A td = td (L) — td{xr) • This lag time 
generally will cause the Mach number to drop before the drive motor and IGV system control inputs can correct for 
it. Thus, to account for this time delay, the drive motor and IGV system control inputs must precede the total pressure 
disturbance in the test section by the same time lag. In effect, this is the essence of the proposed model predictive Mach 
number control strategy whereby the control inputs have to be computed from an accurate mathematical model of the 
wind tunnel in a future time prior to the start of a continuous pitch polar when the error signals are not yet available. In 
contrast with the current Mach number feedback control strategy, the model predictive Mach number optimal control 
is a feedforward, open- loop scheme. The overall drive motor and IGV system model predictive control inputs can now 
be computed as 

Va(t) 

Rr ( t ) 



AF a (t + A td) 

A Rr (t + A td) 


10 


American Institute of Aeronautics and Astronautics 



To demonstrate the effectiveness of the proposed model predictive Mach number optimal control, we select a test 
section Mach number M (xt) — 0.6, a test section total pressure pb (%t) — 2116 psf. We conduct a simulation for a 
test model being pitched continuously from 0° to 15° in 30 sec with a total pressure disturbance as shown in Fig. 7. 

The compressor speed is nominally at 455 rpm and the IGV flap deflection is at 11.5°. The geometry information is 
taken from the NASA Ames 1 1-Foot Transonic Wind Tunnel with a length L = 795 ft and the test section at xt = 470 
ft. The period is computed to be th — 15.6 sec, and the time delay between the test section and the compressor inlet 
is computed to be A td — 3.6 sec. The system matrices are G = [ 25.5415 —684.4960 0 ] , H = 1.6013, and 
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The weighting matrices are selected to be P = 0.001, Q = diag (0, 0.01, 0.01), R = diag (l x 10 6 , 1, 0). 



Fig. 7 - Test Section Total Pressure Disturbance 

The results of the control simulation are presented in Figs. 8 to 13. Fig. 8 is a plot of the feedback IGV field 
voltage control perturbation and Fig. 9 is the plot of the drive motor rotor resistance control perturbation for both 
7 = 0 and 7 = 1. The IGV and the drive motor control perturbation inputs for 7 = 1 are greater in amplitude than 
those for 7 = 0 as expected since the time horizon parameter 7 indicates the degree of a control effort which is the 
largest for an infinite time horizon control. The effect of fluid transport delay in the first 3 sec duration is noted when 
the control inputs are essentially zero. The model predictive control scheme according to Eq. (49) in effect is designed 
to remove this time delay in the control inputs. 

With reference to Fig. 10, the test section Mach number responses for various control schemes are computed from 
the fully nonlinear model. The feedback control schemes for both 7 = 0 and 7=1 clearly are not able to maintain the 
Mach number to within the tolerance band during the dismrbance generated by the continuous pitch of the test model 
in the first 15 sec. Because of the time delay during which the feedback control inputs are inactive as seen in Figs. 8 
and 9, the Mach number drops below the tolerance band before the error signals become available from the sensors at 
the compressor to allow the feedback control inputs to be computed. Once the feedback control inputs are applied to 
the drive motors and the IGV system, the Mach number begins to stabilize and cause the Mach number to eventually 
return to its desired set point 
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Fig. 8 - IGV Field Voltage Perturbation 



Fig. 9 - Drive Motor Rotor Resistance Perturbation 


In contrast, the model predictive control scheme apparently is much more effective in handling the total pressure 
disturbance. The model predictive control scheme for 7 = 0 is an improving control relative to the feedback control 
schemes since it reduces the deviation in the Mach number to within the tolerance band, but suffers a steady state error. 
In contrast, the model predictive control scheme for 7 = 1 is clearly superior to all other schemes in that it is able to 
hold the Mach number to well within the tolerance band throughout the time horizon and converges to the set point 
faster than all other control schemes. This clearly demonstrates the effect of the time horizon parameter 7 . We also 
note that if there were no coirective control, the Mach number would have continued to fall off untilit would'have 
reached a new steady state value. 
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Fig. 10 - Mach Number Response 



Fig. 1 1 - Compressor Speed Response 

Fig. 1 1 shows the compressor speed responses to various control schemes. The model predictive control scheme 
for 7 — 1 provides a much better regulation than all other schemes. Due to the integral action, the compressor 
speed error is virtually right on the set point and well within the desired tolerance of ±0.1 rpm. Both the model 
predictive control scheme and the feedback control scheme for 7 = 0 appear to suffer a.sHghdy_p_oor regulation as the 
compressor speed error does not tend to zero. This is obviously a result of insufficient integral control efforts for a * 
short time horizon control. It is also noted that if there were no corrective control, the compressor speed would have 
continued to drop until a new steady state set point would have been established. 

The IGV flap deflection responses are presented in Fig. 12. The differences between the model predictive control 
scheme and the feedback control scheme illustrate the effect of the time delay quite clearly. Since the compressor speed 
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is not sensitive to various control schemes, it can be seen that the test section Mach number response is predominantly 
controlled by the IGV flap deflection. 



t, sec 


Fig. 12 - IGV flap deflection Response 



t, sec 

Fig. 13 - Mach Number Nonlinear Perturbation Response 

To demonstrate the effectiveness of the model predictive control scheme for 7 = I over the subsonic operating 
envelope, we perform a control simulation for all Mach numbers from 0.4 to 0.9. The results are plotted in Fig. 13. 
As. _can be. seen, the predictive control scheme for 7 = 1 is highly effective for all Mach numbers up to 0.8. The 
Mach number perturbation is well within the tolerance of ±0.001, However, at a Mach number of 0.9, there is a 
rapid increase in the Mach number perturbation that exceeds the tolerance. However, the largest perturbation is about 
—0.0018 which may be acceptable. The sudden change in the Mach number perturbation at a Mach number of 0.9 
is directly attributed to the well-known phenomenon of transonic flow where many linear perturbation theories are 
broken down near a Mach number of unity. 


14 

American Institute of Aeronautics and Astronautics 





It should be noted that in the above model predictive control simulation, we assume a perfect knowledge of the drag 
coefficient of the test model. In reality, this information is generally not necessarily known in advance. Therefore, in 
order to implement a predictive control, a parameter estimation process must be incorporated to provide an estimation 
of the disturbance. In general, the accuracy of a predictive control is predicated upon the goodness of the estimation 
of the disturbance. Any significant deviation from the true value of the estimated parameter of the disturbance will 
likely cause the Mach number to exceed the tolerance band. Since the total pressure disturbance is a function of the 
drag force experienced on the test model, it may be readily estimated from a typical wind tunnel measurement of the 
three aerodynamic force components, namely; axial force coefficient Ca> side force coefficient Cs, and normal force 
coefficient Cjv. These force components can be resolved in a direction parallel to the test section longitudinal axis 
using the Eulerian angles of pitch ip, yaw <p, and roll (p from which the drag coefficient Cd can be estimated according 
to 

Cd — C a cos ip cos <p 4* Cs (cos ip sin <p cos cp — sin ip sin tp) 4 - CW (cos ip sin <p sin <p 4 sin ip cos <p) (49) 

During a pitch polar, the aerodynamic forces generally vary as functions of the pitch angle as well as the Mach 
number and the Reynolds number. A database of the aerodynamic forces can be established as a table lookup process 
for various combinations of the Mach number, the Reynolds number, and the pitch angle. In lieu of the lookup table, 
the aerodynamic force relationships can also be established via an adaptive parameter estimation process such as a 
recursive least-square or a neural network algorithm which provides a functional approximation of the aerodynamic 
force dependency. Using the knowledge of the time response of the model support system, a trajectory of the total 
pressure disturbance parameter w (t) can then be predicted for a given pitch polar. The trajectory of the total pressure 
disturbance is then used to compute the corrective optimal control inputs for the model predictive control scheme. 
Using the model predictive results of the optimal control, the compressor control is switched to an open-loop mode 
and begins its actuation simultaneously with the model support system. At the end of the feedforward mode, the 
compressor control is then switched back to a feedback mode. The model predictive Mach number control scheme 
thus potentially offers a significant advantage over the conventional feedback approach which has been known to be 
unable to regulate the Mach number during a continuous pitch polar in the NASA Ames 1 1-Foot Transonic Wind 
Tunnel 


Conclusion 

This paper describes a new model predictive optimal control method for the NASA Ames 1 1-Foot Transonic Wind 
Tunnel. The wind tunnel is modeled by the 1-D unsteady Euler equations whose boundary conditions prescribe a 
controlling action to the wind tunnel by an empirical compressor performance model. Furthermore, the compressor 
performance in turn is controlled by the drive motors and the inlet guide vane system whose dynamics are prescribed 
by a set of ordinary differential equations. In effect, the Euler equations become coupled to the ordinary differential 
equations via the boundary conditions. This wind tunnel math model has demonstrated an excellent agreement with 
test data. An adjoint method is used to establish necessary conditions of optimality of this type of systems. The 
theoretical results are then used to develop a model predictive optimal control that can be computed by solving the 
Riccati matrix equations. By shifting the control perturbation inputs ahead by a time delay between the test section 
and the compressor inlet as computed by the math model, the model predictive control is able to maintain the test 
section Mach number during a continuous pitch motion as demonstrated by a simulation study. The results of this 
study shows that the model predictive Mach number control potentially offers a much better disturbance rejection than 
a Mach number feedback control. 
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